Astronomy & Astrophysics manuscript no. zanni jd 


© ESO 2008 


February 5, 2008 





MHD simulations of jet acceleration from Keplerian accretion disks 

The effects of disk resistivity 

C. Zanni 1 ' 2 , A. Ferrari 2 ' 3 , R. Rosner 3 , G. Bodo 4 , and S. Massaglia 2 



o 
o 

(N 



Laboratoire d' Astrophysique de Grenoble, 414 Rue de la Piscine, BP 53, F-38041 Grenoble, France 
e-mail: Claudio . Zanni@obs . u j f -grenoble . fr 

2 Dipartimento di Fisica Generale dell'Universita, via Pietro Giuria 1, 10125 Torino, Italy 
e-mail: ferrari@ph.unito.it 

3 Department of Astronomy and Astrophysics, University of Chicago, 5640 S. Ellis Av., Chicago, IL 60637, USA 
e-mail: r-rosner@uchicago.edu 



4 Osservatorio Astronomico di Torino, Viale Osservatorio 20, 10025 Pino Torinese, Italy 
e-mail: bodo@to . astro . it 



in 



Received ; accepted 

ABSTRACT 



Context. Accretion disks and astrophysical jets are used to model many active astrophysical objects, viz., young stars, relativistic 
stars, and active galactic nuclei. However, extant proposals on how these structures may transfer angular momentum and energy from 
, disks to jets through viscous or magnetic torques do not yet provide a full understanding of the physical mechanisms involved. Thus, 

global stationary solutions do not permit an understanding of the stability of these structures; and global numerical simulations that 
include both the disk and jet physics have so far been limited to relatively short time scales and small (and possibly astrophysically 
unlikely) ranges of viscosity and resistivity parameters that may be crucial to define the coupling of the inflow-outflow dynamics. 
, Aims. Along these lines we present in this paper self-consistent time-dependent simulations of supersonic jets launched from mag- 

netized accretion disks, using high resolution numerical techniques. In particular we study the effects of the disk magnetic resistivity, 
parametrized through an a-prescription, in determining the properties of the inflow-outflow system. Moreover we analyze under which 
conditions steady state solutions of the type proposed in the self-similar models of Blandford & Payne can be reached and maintained 
in a self-consistent nonlinear stage. 
Methods. We use the resistive MHD FLASH code with adaptive mesh refinement, allowing us to follow the evolution of the structure 
5^ ' for a time scale long enough to reach steady state. A detailed analysis of the initial configuration state is given. 

Results. We obtain the expected solutions in the axisymmetric (2.5 D) limit. Assuming a magnetic field around equipartition with the 
thermal pressure of the disk, we show how the characteristics of the disk-jet system, as the ejection efficiency and the energetics, are 
affected by the anomalous resistivity acting inside the disk. 
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1. Introduction some basic models have been proposed and shown to be suc- 
cessful (for rather recent reviews see Pudritz et al. 2006a| for 
Astrophysical jets are an important component in many active ySO and Ferrari EH |2004| for AGN). The overall idea is to 
astrophysical objects, from young stellar objects (YSO) to rela- extract energy and angular momen tum from the accreting matter 
tivistic stars and galactic nuclei. In particular, Herbig-Haro (HH) and t0 feed it into a plasma that is acce l e rated in two opposite 
outflows are detected in stellar forming regions around T-Taun directions along the rotation axis of the disk . Lovelace (fl976b 
stars, charactenzerimainly by optical emission line spectra (e.g. and Blandford (fl976b proposed independently that this can be 
Reipurth & Bally [1996D. On a much larger scale, relativistic ra- done by electromagnetic forces. In particular Blandford & Payne 
dio jets are accelerated in the innermost cores of Active Galactic derived a steady state MHD solution for an axisymmetric 
Nuclei (AGN) (e.g. Giovanmni[2004): their emitting component ma gnetocentrifugally driven outflow from a Keplerian disk; for 
is synchrotron relativistic electrons, with a cold proton compo- the outflow to be launched the poloidal magnetic field lines must 
nent or, most likely, a Poynting flux electromagnetic component be incline d less than 60° with respect to the plane of the accre- 
(De Young|2006|). tion d isk p n me ot h er h and, Sa uty & Tsinganos ( [T994j l, Sauty 
Despite of being characterized by extremely different space, e t al. (120021 , Sauty et al. ( 120041 1 have derived meridionally self- 
time and energy scales, it is commonly accepted that all these similar models to study the launching mechanism from the hot 
systems derive their energy from accretion onto a central object corona of the central object. 
(Livio |T999| : the physical origin of these supersonic outflows 

has been related to the dynamical evolution of magnetized ac- T h e magnetocentrifugal mechanism has be en the subject of 

cretion disks around a deep gravitational well. Even if the ac- a series of nu merical studies (Ust yugova et al. [19991 Ou yed & 

celeration and collimation mechanisms of jets are still not clear, Pudrit z [19971 Krasn opolsky et al. \lW9\ Anderson et al. [20041 

Fendt [20061 Pudritz [5006b) based on ideal MHD simulations 

Send offprint requests to: C. Zanni in which the disk is treated as a boundary condition. On one 
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hand they show how a steady solution can be obtained in a few 
dynamical time scales and how the acceleration, collimation and 
stationarity of the outflow depend on the mass loading from the 
disk and on the magnetic field structure. On the other hand the 
back reaction of the outflow on the disk can not be taken into 
account. 

When the structure of the magnetized accretion disk is in- 
cluded self-consistently in the models, a diffusive mechanism 
must be introduced inside the disk to balance the shearing due 
to differential rotation and the inward advection of field lines. 
Moreover also viscous torques, which can transport angular mo- 
mentum radially inside the disk itself, should be, in principle, 
taken into account. 

Fo r instance, Konigl (fl989] l, Wardle & Konigl ( [19931 and Li 
( 1996 ) studied the structure of an ambipolar diffusion dominated 
disk connecting it with a Blandford & Payne solution at the disk 
surface. Ogilvie & Livio (2001) solved the vertical structure of 
a thin (optically thick) magnetized disk taking into account an 
effective a turbulent viscosity and resistivity. 

In a series of papers including turbulent a resistivity 
(Ferreira & Pelletier 119951 Ferreira 119971 ), viscosity (Casse 
& Ferreira 2000a) and entropy generation (Casse & Ferreira 
2000b) the authors calculated radially self-similar stationary so- 
lutions of accretion-ejection structures. Two important features 
emerged clearly from these latter models: first of all it was shown 
that, in order to balance the magnetic and gravitational com- 
pression on the disk itself, the thermal energy must be around 
equipartition with the magnetic energy inside the disk. With 
these conditions, the vertical thermal pressure gradient is the 
only force that can push the mass at the surface of the disk to 
be accelerated in the outflow. Second, the magnetic torque must 
change sign at the disk surface, in order to extract angular mo- 
mentum from the disk and transfer it to the outflow. This con- 
dition determines a strong constraint on the resistive configura- 
tion: the magnetic diffusivity must be rather high (a ~ 1) and 
anisotropic, the diffusion of the poloidal field being smaller than 
the diffusion of the toroidal component. 

The first time-dependent numerical simulations in which the 
structure of the magnetized accretion disks is included (Uchida 
& Shibata fl985l or more recently Kato et al. 120021 ) showed that 
the interaction between a geometrically thin rotating disk and a 
large scale magnetic field that was initially uniform and verti- 
cal originates a transient state in which a strong toroidal field is 
generated that expels matter in the direction perpendicular to the 
disk plane ("sweeping magnetic twist mechanism"). One of the 
limits of these models is that the short simulated time scales and 
the ideal MHD approximation lead to the formation of a tran- 
sient and highly unstable outflow. Moreover the small density 
contrast between the disk and the surrounding corona assumed 
in Uchida & Shibata ( 119851 ) enables the disk to lose its angular 
momentum on a very short time scale. These works have been 
recently updated in Kuwabara et al. (120051 ) who studied the evo- 
lution of a thick magnetized torus assuming a constant resistiv- 
ity throughout the computational domain: even if the presence of 
such a high and uniform resistivity must be justified, the authors 
showed a quasi-stationary outflow launched from the in the inner 
radii of the torus. 

Up to the present, the best effort to produce an accretion- 
ejection structure recurring to time-dependent simulations has 
been performed by Casse & Keppens (2002, 2004) who showed 
how a quasi-stationary jet can be launched from the equipartition 
regions of a resistive accretion disk. On the other hand, their re- 
sistive configuration, isotropic with a = 0.1 is rather different 
from the one predicted by the self-similar steady models (Casse 



& Ferreira 2000a): despite of the use a stretched grid, the res- 
olution of these simulations is rather low and it is likely that 
numerical dissipative effects are quite important. 

In the present paper we present a numerical study of resis- 
tive MHD axisymmetric accretion-ejection structures performed 
with the high resolution code FLASH using Adaptive Mesh 
Refinement. The aim of the paper is to simulate the disk-jet con- 
figuration over long time scales in order to determine the effects 
of different configurations of an a resistivity, by varying its value 
and its degree of anisotropy, in determining the properties of the 
system. Moreover we want to test whether a stationary state cor- 
responding to the Blandford & Payne self-similar solution can 
be reached and maintained. A critical point in the simulation 
is of course the choice of the initial configuration, in particu- 
lar the magnetic configuration and the equilibrium of the disk. 
We have solved for an analytic self-similar equilibrium config- 
uration of the disk including gravitational, centrifugal, thermal 
pressure and Lorentz forces; this solution defines also the ini- 
tial magnetic field. We do not include physical viscosity: as it 
has been shown both in stationary (Casse & Ferreira 2000a]) and 
time-dependent contexts (Meliani et al. 120061 ). for conventional 
values of the disk turbulence {a ~ 0.1 - 1, Prandtl number ~ 1) 
the viscous torque is less efficient than the magnetic in extracting 
angular momentum from the disk. 

The paper is organized as follows. Section 2 is dedicated to 
illustrating the equations and the numerical code used; the initial 
configuration and the boundary conditions are discussed in de- 
tail. In Section 3 the numerical simulations are presented for dif- 
ferent magnetic diffusivities and anisotropies. These parameters 
are linked to the mass accretion rate, required torque and liber- 
ated accretion energy. The acceleration mechanisms at the disk 
boundary and far above the disk are discussed, also referring to 
the established current circuits. The jet ejection efficiency is dis- 
cussed in Section 4 in relation to the accretion rates. In Section 
5 we discuss the angular momentum transport and in Section 6 
the energy budget of the inflow/outflow system is evaluated. In 
Section 7 we test if one of our cases can be characterized as a 
stationary solution. Comments about the effect of Ohmic dissi- 
pation in the disk are given in Section 8. Finally, Section 9 sum- 
marizes our results and puts them in context with other analytical 
and numerical models. 



2. The numerical model 

2.1. MHD equations 

We model the interaction between an accretion disk and the 
magnetic field which threads it within a resistive MHD frame- 
work. The system of equations that we solve numerically con- 
veys therefore the conservation of mass: 



% + V ■ {pu) = , 
ot 



(1) 



where p is the mass density and u is the flow speed; the momen- 
tum equation: 



dpu 
~~dT 



+ V 



B B\ 
puu + IP + — — | / 



BB 



+ p Vd) g = . 



(2) 



where P is the thermal pressure and B is the magnetic field. This 
equation takes into account the action of thermal pressure gradi- 
ents, Lorentz forces and gravity as determined by the potential 
Og = —GM/ Vr 2 + z 2 representative of the gravitational field of 
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a central object of mass M. The evolution of the magnetic field 
is determined by the induction equation (Faraday's law): 



— + VxE = 0. 

dt 



(3) 



where the electric field E is determined by the Ohm's law 
E — -u x B + f\J . Notice that the equations are written in 
a non-dimensional form, hence without \n and po coefficients. 
The electric current J appearing in the Ohm's relation is deter- 
mined by the Ampere's law J = V x B. The magnetic resistivity 
f\ is indicated as a two-tensor to take into account anisotropic 
diffusive effects: in our simulations we will consider a diagonal 
resistivity tensor rjij whose non-zero components are = j] m 
and TJrr = T]zz = 7m- 

Finally the conservation of energy is expressed by 



de 
~dt 



• + P + — I u - (u B) B + fjj x B 



= -A 



,(4) 



where the total energy density 
P 



y- 



i 



pu u B ■ B 

+ — =— + +p® s 



is given by the sum of thermal, kinetic, magnetic and gravita- 
tional energy, y = 5/3 is the polytropic index of the gas. A C00 ] is 
a cooling term defined by the parameter < / < 1 : 



A, 



cool 



*-diss 



(5) 



given by the ratio between the specific radiated energy and the 
Ohmic heating term Ad; ss = f\J ■] . The parameter / therefore de- 
termines the fraction of magnetic energy which is radiated away 
instead of being dissipated locally inside the disk increasing its 
entropy. Finally the system of equations is closed by the equa- 
tion of state of ideal gases P = nKT where n = p/m v (m p being 
the proton mass) is the number density of the gas, T is its tem- 
perature and K is the Boltzmann constant. 

To solve the resistive MHD system of Eqs. ([TJUi we em- 
ploy a modified version of the MHD module provided with 
the public code FLASH3 (Fryxell et al. EOOOi i developed at the 
ASC FLASH Center at the University of Chicago, adopting its 
Adaptive Mesh Refinement (AMR) capabilities. The simulations 
have been carried out in 2.5 dimensions, that is in cylindrical ge- 
ometry in the coordinates r, z assuming axisymmetry around the 
rotation axis of the disk-jet system. The algorithm implemented 
belongs to the class of high resolution Godunov schemes which 
are the best suited to study supersonic flows: we therefore used a 
linear reconstruction of primitive variables with a minmod lim- 
iter on pressure and flow speed and a Van Leer limiter on density 
and on the magnetic field components; second order accuracy in 
time is obtained thanks to an Hancock predictor step on the prim- 
itive variables while, to compute the fluxes needed to update the 
conservative variables, we implemented an HLLE solver, that 
is an approximate linearized Riemann solver, that assumes a- 
priori a two-wave configuration for the solution. To control the 
solenoidality of the magnetic field (V • B — 0) the eight-wave ap- 
proach (Powell et al. 1 19991 1 was used, which is known for simply 
advecting the monopoles, while a parabolic diffusion operator 
(Marder fl987l see also Dedner et al. [2002) was added to the in- 
duction equation to diffuse them. We finally ensured to use an 
angular momentum conserving form of the <p component of the 
momentum equation Eq. (f2]i and a poloidal current flux conserv- 
ing form of the <p component of the induction equation Eq. (01. 



FLASH is freely available at http://flash.uchicago.edu 



2.2. Initial conditions 

In the initial setup of our simulations we model a disk rotat- 
ing with a slightly sub-Keplerian speed threaded by an initially 
purely poloidal magnetic field. The disk initial model is derived 
imposing equilibrium between the forces initially intervening in- 
side the disk, namely, gravity, centrifugal force, thermal pressure 
gradients, and Lorentz force; the disk setup is therefore a solu- 
tion of the following system of equations: 



dP 

dP 
or 



dz 
dr 



— Jd,B r 



+ J*B- + 



(6) 



(7) 



The system of Eqs. (|6]|7]) can be easily solved in separable 
variables assuming radial self-similarity; with this assumption 
all the physical quantities U are given by the product of a power 
law for r with a function of the variable x — z/r 



(8) 



where z — corresponds to the midplane of the disk, fu is 
an even function such that fu(0) = 1 or an odd function such 
that fu(0) = 0, depending on the symmetry of the variable with 
respect to the midplane of the disk. For the physical quantities 
which have an even symmetry (P, p, u r , u$, B z ) Uq therefore 
represents their value at r = ro, z- 0. 

Self-similarity requires that all the characteristic speeds, 
namely sound and Alfven speed, and flow speeds should scale 
as the Keplerian velocity (oc r~'/ 2 ) on the midplane of the disk. 
Imposing also a polytropic relation between the disk density and 
pressure, that is P — Pq (p/po) y , the power law coefficients f3u 
are determined as follows: 



B»t=B Ur =B Ui = -1/2 
B Br = Pb z = -5/4 



Bp = -5/2 
B P = -3/2 . 



The initial poloidal magnetic field has been set through the 
flux function *P to ensure the solenoidality of the field: 



3/4 



,5/4 



4 2 / r \ m " 



(9) 



The components of the field, obtained thanks to the simple rela- 
tions: 

1W 1 <9¥ 
B z = -— B r = — 

r or r oz 

obviously fulfill the self-similarity requirements. The parameter 
m determines the height scale on which the initial magnetic field 
bends, where a value m — > oo gives a perfectly vertical (B r = 0) 
field. 

Given the poloidal magnetic field, the vertical equilibrium 
Eq. (|6]l can be numerically solved to give the vertical profiles 
of disk density (f p ) and pressure (fp = fl). On the other hand 
the radial equilibrium Eq. (0 can be solved to determine the 
disk rotation speed u^. The solution will therefore depend on the 
following non-dimensional parameters: 



p r 
pGM 



z=0 



(10) 



jP/p and the 



which is the ratio between the sound speed c s 
Keplerian rotation speed Vk = y/GM/r evaluated on the disk 
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midplane: this quantity determines the disk thermal height scale 
H through the relation H — er (see Frank et al. 120021 Ferreira & 
Pelletier [l995l ); the magnetization parameter: 



B 2 



(11) 



which gives the ratio between the magnetic and thermal pressure 
evaluated on the disk midplane. The disk rotation velocity at the 
midplane is slightly sub-Keplerian, with the deviation depending 
on the parameters (e, /u and m): 



U <P\~=> 



=o 



2e> I - + „ , 

H4 3m 2 



GM 
r 



(12) 



The toroidal velocity decreases towards the disk surface and it 
falls down to zero where the radial Lorentz force balances the 
gravitational pull and the thermal pressure gradient. 

For the components of the magnetic diffusivity tensor f\ act- 
ing in the disk we adopt an a prescription (Shakura & Sunyaev 
119731 1 in the same vein of Ferreira (1997 and related works) and 
Casse & Keppens (2002 2004). Despite of being a common way 
to parametrize the transport coefficients inside an accretion disk, 
the origin of this anomalous diffusivity is still debated. One of 
the most promising hypothesis states that the anomalous trans- 
port is of turbulent origin, triggered by some disk instability, be- 
ing the magneto-rotational instability (MRI, Balbus & Hawley 
1998) the most accredited one. The tj^ = rj m component is 
parametrized as follows: 



;/„, = a m V A | z=0 #exp|-2^2 



(13) 



where a m is a constant parameter, VaI^o 
the Alfven speed calculated on the disk midplane and H = 
(c s /Qk)I z =o is the thermal height scale of the disk. In the simula- 
tions both the Alfven speed and H are allowed to evolve in time. 
The other components of the diffusivity tensor r\ rr = rj zz = rj' m 
are assumed to be proportional to rj m through an anisotropy pa- 
rameter Xm, which is the inverse of the analogous parameter in- 
troduced by Ferreira & Pelletier (1995): 



(14) 



A ratio Xm — 1 indicates an isotropic resistive configuration. 
We recall that the presence of an effective resistivity inside the 
disk allows the magnetic field to break the "frozen-in" condition 
and the matter to slip through the field lines. The component 
r] m , therefore indicated as poloidal resistivity, allows the flow to 
slip through the field in the poloidal plane while rj' m , indicated as 
toroidal resistivity, controls the diffusion of the toroidal compo- 
nent of the field. 

In a steady situation the accretion motion should be consis- 
tent with the poloidal magnetic field configuration as stated by 
the poloidal induction equation in its stationary form: 



u z B r - u r B z - rj m J$ 



(15) 



which expresses the balance between the advection and the dif- 
fusion of the poloidal magnetic field. Therefore we initially im- 
pose an accretion flow inside the disk solving Eq. (TT3T > with the 
condition u z = ^u r . Due to the self-similarity requirements, the 



radial accretion speed scales initially as the Keplerian speed on 
the midplane of the disk: 



GM 
r 



(16) 



As the poloidal resistivity rj m exponentially decreases towards 
the disk surface, the initial accretion flow cancels out outside the 
disk. 

It must be pointed out that in the initial conditions, since 
there is no toroidal magnetic field, there is no mechanism of an- 
gular momentum transport which can support the initial accre- 
tion flow: the angular momentum transport associated with the 
toroidal field will be triggered by torsional Alfven waves due 
to the differential rotation between the midplane and the surface 
of the disk. Moreover, since we will assume a magnetic field 
around equipartition with the thermal energy on the midplane 
of the disk (see Section [231 ), the time scale on which the trans- 
port of angular momentum will become effective, given by the 
Alfven crossing time of the disk thickness, is comparable to the 
local period of rotation of the Keplerian disk and to its epicyclic 
frequency. 

On top of the disk we prescribe an hydrostatic spherically 
symmetric atmosphere for which we impose, as for the disk, a 
polytropic relation between pressure P a and density p a : P a = 
^ao (Pa/Pao) r - The density and pressure distribution will be there- 
fore given by: 



Pa = PaO 



ffl 



■ PaO" 



7 — 1 GM I ro \ r-i 



7 ro 



(17) 



where p a o is the value of the atmosphere density at the spherical 
radius R = ro- The initial position of the disk surface is located 
whereas the disk and atmosphere pressures are equal. The initial 
position of the disk surface is therefore determined by its tem- 
perature, defined by e, by the field intensity fi, its inclination, 
given by the parameter m and by the density contrast between 
the disk and the corona p a o/po- 



2.3. Units and normalization 

Since in the formulation of the problem we did not introduce 
any specific physical scale, the system of Eqs. (Q]|4| and the ini- 
tial conditions presented in Section 12.21 can be normalized in 
arbitrary units: the results will be therefore presented in non- 
dimensional units. 

Lengths will be given in units of ro, corresponding approx- 
imatively to the inner truncation radius of the disk; speeds will 
be expressed in units of the Keplerian speed Vko = y/GM/ro 
at r = ro and the densities in units of po, that is the disk initial 
density at r — ro, z — 0. Assuming for ro the following units, 
appropriated for YSO or AGN systems, 



r = 0.1 AU (YSO) 

= lOflschw = KT^ioty pc (AGN) 



(18) 



where /?schw = IGMjc is the Schwarzschild radius, the 
Keplerian speed Vko is given by 



VK0 = 94(£) 1/2 (^Pkms 



1/2 



6 - 7xi ° 4 (K^r I/2kms 



(YSO) 
(AGN) 



(19) 
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Correspondingly time will be expressed in units of fo = ro/VKO- 
* = 1 - 7 (^r' /2 (a^u) 3/2 d ^ (YSO) 



= ^{^{wt^f 2 ^ (AGN) 



(20) 



Expressed in units of fo the Keplerian period of rotation at the 
inner radius of the disk ro is equal to 2n. 

Finally the normalization density po can be chosen by deter- 
mining a suitable mass accretion rate unit Mo = ^qPoVko: 



^ = 3xl0-(-^)(a' /2 (^) 3/2 M yr- 
= 9 (io-i2gcm- 3 )(TOTo) (lofiw) M oyr _I 



(21) 



where, as customary, the first expression corresponds to YSO 
and the second to AGN objects. The values of accretion and out- 
flow mass rates will be presented in non-dimensional units and 
they must be multiplied by the Mo factor of Eq. ( |2"TT i to obtain 
their physical value. Similarly torques and powers which will be 
shown subsequently will be given in units of Jq - r^poV^D an< ^ 
Eo = r oA>^KO respectively: 



Jo- 3 x 



= 1 - 2 >< 1o51 (ich^)(ic^) 3 (to^) 2 dynecn 
*o = 1-9 X 1033 eig ^ 

' x 1046 ( Hi-' 2 gem- 3 ) ( To% ) (loiw) 



(22) 



(23) 



■ 2.6: 



erg s 



Again, the first equations must be used for YSO while the second 
for AGN systems. 

2.4. Boundary conditions 

Besides of restricting our study to axisymmetric structures we 
also assume that the system is symmetric with respect to the mid- 
plane of the disk. The computational domain therefore covers a 
rectangular region with a radial extent [0, 40ro] and a size along 
the z direction equal to [0, 120ro]: axisymmetry is imposed on 
the rotation axis r — while planar symmetry is imposed on 
the disk midplane z = 0. Besides of choosing suitable boundary 
conditions on the outer sides of the domain, an inner boundary 
must be placed inside the computational box in order to avoid 
the singularity of the potential well and of the initial setup at the 
origin of the axis: we therefore define a rectangular box with a 
size r x z = [0, ro] x [0, 0.5ro] which is excluded from the com- 
putation and on whose sides boundary conditions must be set. 

On the outer right boundary (r = 40ro) "outflow" conditions, 
that is zero-gradient, are imposed on thermal pressure, density 
and poloidal velocity components. On one disk thermal height 
scale (z < Hr) the value of the radial velocity in the ghost zones 
is set to the extrapolated value only if negative and zero other- 
wise: in this way we avoid any outflow of matter without im- 
posing the mass accretion rate, which will be determined by the 
dynamical evolution of the system. For u^, B z , the continuity 
of the first derivative is required too. The condition on B r is de- 
termined by imposing the solenoidality of the field, V • B = 0, 
on the last cell of the domain with a first order approximation. 



On the outer upper boundary (z — 120ro) "outflow" condi- 
tions are imposed on all the variables except for B : whose bound- 
ary values are determined to satisfy the V • B — requirement. 

We decided to prescribe the continuity of the first deriva- 
tive of some variables in the radial direction since the power- 
law behavior of the initial conditions along r (see Eq. ([8])) de- 
termines strong gradients of the physical quantities on the right- 
most boundary. This boundary condition affects above all the 
radial component of the Lorentz force F r : 



F r = - 



l J3 

2 dr 



1 dB\ 8B r 

2 dr dz 



(24) 



and therefore the collimation of the outflow. First of all it is easy 
to see that a zero-gradient condition on the toroidal magnetic 
field component B^ cancels the pressure gradient of the magnetic 
pressure (firs term in the right hand side of Eq. (l24l ) while the 
pinching force —B % Jr (second one) is still present. An outflow 
condition on B^ on the right outer boundary thus enforces the 
collimation due to the toroidal field. This issue has been widely 
discussed by Ustyugova et al. (1999) who proposed to use a 
"force-free" condition, where the poloidal electrical current is 
parallel to the magnetic field, in order to remove artificial forces 
originating on the boundaries. Seen from the point of view of 
electrical currents, a zero-gradient condition on Z?^ corresponds 
to have a negative collimating current component J z flowing 
along the right boundary, while stationary models of accretion- 
ejection structures show that on the outer part of the disk the 
current outflows from it thus pushing the flow along the field 
lines without collimating it (see Fig. 13 in Ferreira [1997). The 
continuity of the first derivative of Z?^ on the rightmost bound- 
ary allows to develop a gradient of magnetic pressure associated 
with Z?^ which can counteract the toroidal pinch and generate an 
outflowing positive current J z . 

Even if the effects are less pronounced, also the continuity 
of the first derivative of B z on the rightmost boundary affects 
the collimation of the structure: a pressure gradient directed out- 
wards associated with B z (third term in the right hand side of 
Eq. (T24l ) counteracts the poloidal field tension (fourth term) thus 
producing slightly less collimated structures. 

On the other hand, we noticed that the choice of an "out- 
flow" condition for the toroidal field on the upper boundary did 
not affect the behavior of the outflow as much as the condition 
on the right boundary. We did not notice a huge difference be- 
tween the "outflow" condition, where the poloidal current has 
only a z component, and the "force-free" condition proposed by 
Ustyugova et al. (1 19991 1, where the poloidal current is parallel to 
the poloidal field. 

On the inner boundaries located on the edges of the rect- 
angular region r x z — [0, ro] x [0,0. 5ro], we adopted a simi- 
lar strategy: on the r = ro side we extrapolated all the physical 
quantities imposing the continuity of the first derivative except 
for the poloidal components of the velocity field, for which a 
zero-gradient condition was used, and for the radial component 
of the magnetic field, which is required to fulfill the V ■ B — 
condition; on the z = 0.5ro side we imposed an "outflow" condi- 
tion on all the variables except for the z component of the field, 
determined by imposing the solenoidality of the field. 

All the details on the simulation performed are given in the 
next Section. 
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Table 1. Initial parameters of the simulations performed: poloidal resistivity parameter a m , anisotropy of magnetic resistivity Xm, 
cooling parameter /, initial accretion rate M/Mo, torque J/Jo needed to support the accretion rate between r; = ro and r e = 10ro, 
accretion energy E/Eq liberated between the same radii, equivalent resolution of the adaptive grid. 



Simulation 


a m 


Xm 


/ 


M/M 




J/Jo 


E/E 




Resolution 


1 


0.1 


l 


1 


7.3 x 10' 


■■} 


1.4 x 10~ 2 


3.3 x 10" 


-3 


512 x 1536 


2 


0.1 


l 


1 


7.3 x 10 


-3 


1.4 x 10~ 2 


3.3 x 10" 


-3 


128 x 384 


3 


0.1 


l 





7.3 x 10" 


-3 


1.4 x 10~ 2 


3.3 x 10" 


-3 


512 x 1536 


4 


1 


l 


1 


7.3 x 10" 


-2 


1.4 x 10-' 


3.3 x 10" 


-2 


512 x 1536 


5 


1 


3 


1 


7.3 x 10" 


-2 


1.4 x 10~' 


3.3 x 10" 


-2 


512 x 1536 


6 


1 


3 





7.3 x 10" 


-2 


1.4 x 10-' 


3.3 x 10 


-2 


512 x 1536 



2.5. The simulations 

Once they are normalized with the units given in Section 12.31 
the initial conditions presented in Section l2~2l depend on 6 non- 
dimensional parameters: the ratio between the sound speed and 
the Keplerian speed at the disk midplane e; the disk magneti- 
zation parameter /i; the magnetic height scale of the initial field 
m; the strength of magnetic diffusivity a m and its anisotropy co- 
efficient Xm\ the ratio between the initial atmosphere and disk 
densities p a o/po- All the free parameters except those describing 
the diffusive properties of the disk will be the same for all the 
simulations. We therefore assume a parameter e = 0.1, which 
fixes the initial thermal height scale and temperature T on the 
midplane of the disk, that is: 

T\ z=0 = ^ = ^(^(M-)^)- 1 K 

(25) 

= 5xl0 9 (^) 2 ( T?i ^)- 1 K 

The initial structure of the magnetic field is determined by fix- 
ing the magnetization parameters fi - 0.3 and m = 0.35, which 
gives therefore a magnetic height scale 3.5 times higher than the 
initial thermal height scale of the disk. We recall that a value of 
magnetization around equipartition or slightly below it is gener- 
ally required both in numerical (see Zanni et al. 120041 Casse & 
Keppens [2002 ) and analytical (Ferreira |1997l Casse & Ferreira 
2000a ) modeling of accretion disks launching jets; this condition 
is determined by the equilibrium between thermal pressure gra- 
dients, which vertically support the disk and are responsible for 
the initial mass loading the field lines, and Lorentz forces, which 
tend to pinch the disk. The magnetization parameter yu fixes the 
initial value of the poloidal magnetic field on the disk midplane: 

B\ Z=Q = JT&Fp = 

- 2 6 ( JL\ m ( p° , V /2 (M.\ V1 r^_r 5/4 G 
^■"10.37 vo.i / ^ io 12 g cm-* ] \m ) Io.iau; (26) 

= 1.8xl03(^f(-)( TF ^) 1/2 ( T ^)- 5/4 G 

Finally we assumed a density ratio between the corona and the 

diskp a o/po = 10~ 4 . 

In order to investigate the effects of magnetic resistivity we 
performed a series of simulations varying the value of the a m 
parameter and the anisotropy factor Xm- The summary of all the 
simulations performed is given in Table Q] 

Isotropic resistive configurations (x m = 1) have been studied 
for two different values of the a m parameter: a m = 0.1 (simu- 
lation 1), which corresponds to the magnetic resistivity adopted 
by Casse & Keppens (2002, 2004), and a m — 1 (simulation 4). 
To determine the effects of an anisotropic resistivity, required by 
the steady models of Ferreira & Pelletier (1995), we performed 



a simulation characterized by a m = 1 and^ m = 3 (simulation 5): 
the parameters a m , x m , ^ an d e of this simulation are typical of 
cold self-similar solutions found by Casse & Ferreira (2000a). 

For these standard simulations the adaptive mesh provided 
with FLASH is set up with 7 levels of refinement based on blocks 
of 8x8 square cells, giving an equivalent resolution of 512x1536 
points: this resolution is kept fixed in the disk region (z < 3er) 
and around the inner boundary rectangle, while the grid is free to 
change and adjust the resolution in the outflow region. In order 
to determine the importance of numerical diffusive effects in our 
simulations, we repeated the case characterized by an isotropic 
a m = 0.1 with a resolution four times lower than the usual one, 
thus increasing the numerical dissipative effects (simulation 2): 
for this simulation the adaptive grid is allowed to reach a maxi- 
mum equivalent resolution of 128 x 384 points, while maintain- 
ing the higher one (512 x 1536 points) around the central inner 
boundary. The resolution of this test case is similar to the one 
adopted by Casse & Keppens d2002ll2004b . 

This first set of four simulations is characterized by a cool- 
ing factor / = 1 (Eq. (O) which implies that all the magnetic 
energy Ohmically dissipated is radiated away. We therefore re- 
peated simulations 1 and 5 assuming a cooling factor / = to 
study the effects of the Ohmic heating in the case of a lower 
(simulation 3) and higher resistivity (simulation 6). The results 
presented in Section 07] will refer to the "cold" cases charac- 
terized by / = 1 while a few comments on the "heated" / = 
simulations will be done in Section[8] 

In the fifth column of the Table we show also the mass accre- 
tion rate that we imposed at the beginning of our simulations: we 
recall that this initial rate is determined by balancing diffusion 
and advection of the poloidal field inside the disk. Obviously en- 
ergy and angular momentum must be extracted from the accre- 
tion flow to support this accretion rate: in the sixth and seventh 
columns of the same Table we show also the torque j and the 
power E liberated in the accretion between r; = ro and r e = 10ro. 
The values of j and E are with a good approximation given by: 

j = M(y/GMr e - y/GMn) (27) 
. IGM GM\ 

E - M [^--^J (28) 

3. Production of jets from magnetized accretion 
disks 

All the simulations were carried on up to a time t = 400 which 
corresponds to ~ 63 periods of rotation of the disk at its inner 
radius. On the other hand the final time t = 400 corresponds 
to only 0.25 rotations at the outer boundary r = 40ro. It is very 
unlikely that the outer part of the disk has reached an equilibrium 
stage at the end of the simulation. Therefore the study of the 
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t = 0. t = 50. t = 100. t = 200. t = 400. 




r r r r r 



Fig. 1. Time evolution of density maps in logarithmic scale of the simulation characterized by (a m = 1, Xm = 3, / = 1). Time is 
given in units of to (see text). In these units the Keplerian period at the inner radius of the disk r = ro is equal to 2n. Superimposed 
are sample magnetic field lines: the distance between the field lines is proportional to the intensity of the field. In the last panel 
(t = 400) are also plotted the critical Alfven (dashed line) and fast-magnetosonic (dotted line) surfaces. 




Fig. 2. Same as Fig.Q~]but for the case characterized by (a m = 0.1, Xm — 1, / = !)■ 



accretion-ejection system will be restricted mainly to the inner 
part of the disk (r < 10r , see Section |4]i and to the outflow 
coming from it. 

In all the cases studied we observed a robust outflow to 
emerge from the underlying accretion disk: the solutions show 
a hollow jet, where the central hole corresponds to the "sink" 
region r < ro. The outflows are not completely collimated at the 
end of the runs, that is some matter is outflowing from the com- 
putational box from the outer cylinder at r — 40ro. Nevertheless 
in all the solutions found the part of the outflow coming from the 
inner part of the disk crosses the Alfvenic and fast-magnetosonic 



critical surfaces inside the domain (see Fig. [T]and|2). Since no 
disturbance produced in the super-fast region of the outflow can 
propagate upstream towards the accretion disk, this condition 
obviously ensures that the outer boundary conditions do not af- 
fect the launching region of this part of the outflow. On the other 
hand perturbations can still propagate in a direction transverse 
to the magnetic field lines: if the fast-Mach cones at the bound- 
aries intersects the computational box, the boundary conditions 
can still affect the radial structure and therefore the collimation 
of the outflow (see Ustyugova et al. 119991 ), as already discussed 
in Section |2~4l 
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To overcome the problem of the influence of boundary con- 
ditions on the launching phase, Krasnopolsky et al. (2003) and 
Anderson et al. (2004) restricted the launching of the wind 
to a narrow region of the inner accretion disk and chose an 
initial magnetic configuration as to contain the entire fast- 
magnetosonic critical surface inside the computational domain. 
The closed shape of the critical surfaces observed in the cited 
papers characterizes those solutions more as "X-winds" (Shu et 
al. |19941 l, that is wide-angle outflows from the corotation radius 
of a stellar magnetosphere, while the almost conical shape of 
the critical surfaces of our solutions is typical of extended disk- 
winds. Moreover, if on one hand it cannot be taken for granted 
that an underlying accretion disk can support the magnetic struc- 
ture and provide the mass load imposed as a boundary condition 
by the authors, on the other hand it has been shown (Ferreira et 
al. 12006 ) that "X-winds" have kinematical properties which are 
inconsistent with observations of T-Tauri microjets. 

The structure of the magnetic field allows to distinguish two 
classes of solutions: the cases characterized by a m — 1 show at 
the end of the computation an "ordered" magnetic configuration 
(Fig. [TJ, while the poloidal field lines are strongly warped and 
distorted if ff m = 0.1 (Fig. [2]). An exception is represented by 
the case performed with a lower resolution: despite of having 
a small resistivity parameter a m = 0.1, it does not present the 
characteristic field inversions (Fig. [3). This anomalous behavior 
can be reasonably ascribed to the higher numerical dissipation 
determined by the lower resolution. 

Other differences can be noticed in Fig. [TJ2] in the more 
diffusive case the jet is asymptotically less dense than in the 
less diffusive one; both outflows become super- Alfvenic (dashed 
line) and super-fast-magnetosonic (dotted line) but in the a m 
0. 1 simulation both characteristic surfaces lie closer to the disk. 

Despite of these morphological differences between the 
more and less dissipative cases, the mechanism which drives the 
outflows from the disk is qualitatively the same: in the follow- 
ing Section we will make a few general considerations about the 
forces which determine both the accretion and the ejection flow. 

3.1. Acceleration mechanism 

Since the accretion-ejection mechanism is mainly magnetically 
driven it is worth writing explicitly the Lorentz forces which are 
acting on the disk-jet system. Following Ferreira ( I1997l l we de- 
compose the Lorentz force F — J x B in the directions parallel 
and perpendicular to the poloidal field: 

*ii = ~M rB *) (29) 

F± = -—V x (rBt) + JfB p 

where B p is the poloidal magnetic field while Vy and V x indicate 
the derivatives parallel and perpendicular to the poloidal field 
respectively. These expressions show clearly that the poloidal 
component of the Lorentz force associated with the toroidal 
field is perpendicular to the isosurfaces rB^ = const, which 
are the surfaces along which the poloidal electric current flows. 
Moreover they show that the component parallel to the poloidal 
field F\\ and the toroidal on e F^ ar e linked by the simple relation 
(see also Casse & Keppens |2002ll2004l ): 

*ll = ~F* (30) 




-20 20 



Fig. 3. Density map (logarithmic scale) at t — 400 of the same 
simulation shown in Fig. [2] performed with a resolution four 
times smaller. 

(notice that the toroidal field Z?^ assumes negative values in our 
simulations). This clearly shows that a poloidal current-field 
configuration that accelerates the outflow along the field lines 
(Fy > 0) is accelerating the plasma also in the toroidal direction 
(F$ > 0) thus providing an additional centrifugal force. This 
represents the core of the so called magneto-centrifugal mecha- 
nism: the magnetic energy stored in the toroidal field at the base 
of the outflow both accelerates the plasma along the field lines 
and increases its angular momentum thus providing a centrifugal 
acceleration. As stated by the above relation the relative impor- 
tance of the two mechanisms, poloidal and centrifugal acceler- 
ation, is given by the ratio LbJ /B p : when the toroidal field is 
stronger than the poloidal one, the gradient of B<f, along the field 
lines is the main accelerating mechanism while if IfiJ /B p < 1 
the plasma tends to corotate with the mainly poloidal magnetic 
field and the centrifugal force is dominant. 

Inside the accretion disk the conditions are reversed: as in 
a unipolar inductor, a radial electric current develops inside the 
conducting plasma, thus providing a toroidal force -J r B, which 
acts to slow down the rotation. According to Eq. (l30l , since 
the Lorentz force brakes the matter in the toroidal direction, a 
poloidal pinching also occurs. Therefore, since the disk is ver- 
tically pinched by the gravity as well as by the magnetic pres- 
sure, only the thermal pressure gradient can ensure a quasi-static 
vertical equilibrium gently lifting the accreting matter towards 
the surface of the disk where it can be accelerated to form the 
outflow. As pointed out in Ferreira ( 119971 1, the transition from 
accretion to ejection can be therefore achieved if the sign of the 
magnetic torque F^, changes sign at the disk surface in order to 
provide a magnetic acceleration. 

This scenario is confirmed by our simulations: in Fig.|4]we 
plot, at the end of our simulations, the total force, given by the 
sum of gravity, Lorentz force and thermal pressure gradient, act- 
ing along the z direction on the disk height scale; these curves 
are obtained by averaging inside the region ro < r < 10ro for the 
cases (a m = 1, Xm — 3, / = 1, solid line), (a m = I, Xm = L 
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z/H 

Fig. 4. Total force acting on the disk height scale along the z 
direction for the cases (a m = 1, Xm = 3, / = 1, solid line), 
(a m — 1 , Xm = 1, / = 1 . dotted line) and (a m — 0. 1 , Xm — 1, / = 
1, dashed line). Along the curves are also indicated the points 
where the magnetic torque (triangles) and the z component of the 
Lorentz force(squares) change sign. The curves are calculated at 
f = 400. 
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Fig. 5. Average ratio between the toroidal and poloidal mag- 
netic fields along the outflows of the simulations characterized 
by (a m = l,Xm = 3, / = 1, solid line), (a m = l,Xm = h f = 1, 
dotted line) and (a m = 0.1, Xm — 1, / = 1, dashed line). The 
curves are calculated at t = 400. 



/ = 1, dotted line) and (a m — 0.1, Xm — 1> / = 1> dashed 
line). Along each line we also indicated the location of the points 
where the magnetic torque (triangles) and the total z-component 
of the Lorentz force (squares) change their sign. It is possible to 
notice that in all the solutions shown the total force changes its 
sign when the disk is still pinched and braked by the magnetic 
field, below both the triangles and the squares: it is therefore the 
vertical thermal pressure gradient that provides the first vertical 
acceleration turning the accretion motion into an outflow; the 
magnetocentrifugal mechanism becomes effective only in corre- 
spondence of the triangles, where the plasma is accelerated both 
in the toroidal direction and along the poloidal field lines. 



Important differences can be noticed between the three cases 
shown: decreasing the value of the poloidal or the toroidal mag- 
netic diffusivity the location of the points where the total force 
and the magnetic torque change their sign are located at lower 
height scales. A possible explanation of this behavior can be 
found in Ferreira (1997): to change the sign of the magnetic 
torque it is in fact necessary that the radial current J r responsible 
for the braking of the disk decreases vertically on a disk scale 
height; it has been shown that in a stationary situation the radial 
current falls off on a height scale proportional to V^mAfm (see 
Eq. (Bl) in Ferreira fl 9971 >. Consistently with our simulations, 
the change of sign of the torque therefore occurs at higher height 
scales increasing the value of a m or the anisotropy factor Xm- A 
second noticeable difference is that while in the two higher dif- 
fusivity cases (a m — 1) the torque changes sign when the disk is 
still pinched by the poloidal magnetic pressure (squares), in the 
a m = 0.1 simulation it changes above this point: in this case the 
z-component of the Lorentz force provides an additional source 
of mass loading before the magnetocentrifugal effect becomes 
effective. 

Another important difference between these three simula- 
tions is shown in Fig. [5] where we plotted a radial average of the 
ratio between the toroidal and poloidal fields \BJ/B P along the 
outflows of the three cases: according to Eq. ( 1301 . these curves 
show that the centrifugal effects are stronger when the magnetic 
diffusivity is anisotropic (solid line), while for a low diffusivity 
(dashed line) the toroidal field pressure gradient is the domi- 
nant accelerating force. Moreover the increase of this ratio with 
z shows that while near the disk surface the centrifugal force is 
effectively contributing to the acceleration of the outflow, it be- 
comes negligible far from it. 



3.2. Current circuits 

A suitable way to understand how the accretion-ejection mech- 
anism works, is to analyze the circulation of the poloidal cur- 
rent in the disk-outflow system: in Fig. [6] we show the poloidal 
current circuits (solid thin lines) superimposed to density maps 
in logarithmic scale at t = 400. Sample field lines (solid thick 
lines) and the poloidal speed vectors are also plotted; the left 
panel refers to the (a m = 0.1, Xm = 1, / = 1) case while the 
right one to the (a m - l,Xm = 3, / = 1) simulation. The current 
circuits are given by the rB^ = const, isosurfaces: the poloidal 
currents circulate counter-clockwise and the Lorentz forces are 
directed outwards perpendicularly to the closed circuits. As dis- 
cussed before, a strong radial positive current flows along the 
disk midplane, thus providing a braking and vertically pinching 
force; at the disk surface the current changes direction until a 
force component accelerating the plasma along the poloidal field 
appears. 

Some differences can be noticed between the two panels: 
due to a strong advection determined by the accretion flow, the 
poloidal field lines are much more inclined inside the disk in the 
less diffusive case (left panel) thus providing a vertically directed 
magnetic tension which mass loads the outflow as it was shown 
in Fig.|4](see the position of the square symbol along the dashed 
line); Moreover, it is possible to see that the different shape of 
the current circuits, together with the higher inclination of the 
poloidal field lines, develops in the a m = 0. 1 a strong vertical 
Lorentz force which strongly bends the field lines, as it was al- 
ready shown in Fig. [2] 

This effect is even more evident if we look at a three- 
dimensional rendering of the two simulations discussed here: in 
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Fig. 6. Poloidal current circuits (solid 
thin red lines) at t — 400 of the cases 
(a m = 0.1, Xm = 1, / = 1, left 
panel) and (a m = 1, Xm = 3, / = 
1, right line). Plotted are also sample 
poloidal field lines (solid thick yellow 
lines) and the poloidal speed vectors. In 
the background density maps in loga- 
rithmic grayscale are shown. 



Fig. |7] we show a 3D rendering of density maps in logarithmic 
scale for the two cases (a m = 0.1, Xm = 1 , / = 1, left panel) and 
(a m = l,Xm = 3,/ = 1, right panel) at a time t = 130. It is pos- 
sible to see that in the more diffusive case the field lines (solid 
blue lines) are gently wrapped around the magnetic surfaces. On 
the other hand, in the a m = 0.1 simulation the footpoints of the 
field lines are advected towards the central object, not balanced 
by a strong enough diffusion; a strong differential rotation along 
the field line therefore appears, the footpoint of the line rotating 
much faster than its opposite end; a strong toroidal field develops 
at the footpoint which gives the strong force directed vertically 
shown in Fig.|6]that completely distorts the field lines. 

The mainly toroidal feature visible in the left panel of Fig. [7] 
propagates vertically along the outflow axis similarly to a "mag- 
netic tower" (Li et al. 120011 Lynden-Bell 2003). Moreover this 
mechanism is not transient, repeating itself on every field line 
when its footpoint is advected towards the faster rotating central 
part of the disk. Since the a m = 0. 1 simulation at a low resolu- 
tion shows a magnetic configuration similar to the right panel of 
Fig. |7l it is possible to speculate that the higher numerical dif- 
fusion can effectively balance the advection of the footpoints so 
that no strong differential rotation along the field lines is present. 

Thanks to the current circuits shown in Fig. [6] we can also 
better understand how the toroidal field can act to collimate the 
outflow against the push of the centrifugal force of the rotating 
jet. It is in fact usually assumed that the self-generated toroidal 
field can automatically ensure the collimation of the outflow 
through the so called "hoop stress" due to the magnetic tension 
of the field and equal to -B^/r. What emerges clearly from both 
the panels of Fig. |6]is that in the outer part of the outflow, where 
the currents flow out from the disk with a positive J z compo- 
nent, the Lorentz force associated with the toroidal field pushes 
the plasma outwards, thus uncollimating the jet: in fact, in this 



situation the outward directed pressure gradient of the toroidal 
field is stronger than the collimating tension. It must be pointed 
out that the shape of the electric current circuits opposing the 
collimation of the outer layers of the outflow strongly depends 
on the outer boundary conditions on B^, as discussed already in 
Section 12.41 On the other hand this behavior is not just a nu- 
merical artifact: the current circuits of our simulations have the 
typical butterfly-like shape of analytical models of disk-winds 
(see Fig. 13 in Ferreira Tl 9971 1 . The outflow starts to be colli - 
mated where the circuit closes back with a negative J z compo- 
nent so that the "hoop stress" exceeds the pressure gradient. As 
the poloidal velocity vectors in Fig.|6]show, in our solutions only 
the inner part of the outflow is almost cylindrically collimated 
inside the computational domain. 



4. Accretion rates and ejection efficiency 

After having analyzed the acceleration and collimation mecha- 
nisms we now take into account the accretion and outflow rates 
which characterize our solutions. In Fig. [8] we plot for all the 
four simulations in which we suppressed the Ohmic heating, the 
accretion rates M a at t — 400, defined as: 

XI. 6H 
pu T dz 
1.6H 

as a function of the cylindrical radius r. H(r) is the thermal height 
scale of the disk defined as H = (c s /Ok)I z =o : the above inte- 
gral is calculated up to 1 .6H, where the magnetic diffusivity de- 
fined by Eq. ( fl~3"l ) becomes negligible and the plasma, now in an 
ideal MHD regime, should flow out from the disk almost par- 
allel to the poloidal magnetic field. It is possible to see that for 
all the curves shown, the accretion rate has a sudden decrease 
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Fig. 7. Three dimensional rendering of density at t — 130 for the cases (a m = 0.1, Xm = 1, / = 1, left panel) and (a m = l,Xm = 3, 
/ = 1, right panel). We also plotted four sample magnetic field lines wrapped around the same magnetic surface. 



for r > 10: this is due to the fact that the outflow is extract- 
ing efficiently mass and angular momentum only for r < 10, 
which is therefore identified as the "launching region". As al- 
ready pointed out in Section [3] at the end of our simulations the 
outer part of the disk has still not reached a dynamically relevant 
time scale and therefore had no possibility to fully develop an 
outflow. Moreover the flow emerging from the inner radii of the 
accretion disk becomes super-fast-magnetosonic inside the com- 
putational domain and therefore the "launching region" will not 
be affected by the outer boundary conditions. We can therefore 
define a control volume delimiting this part of the disk: this vol- 
ume is delimited by an inner cylinder of radius r; = 1 and height 
-1.6H(>i) < z < 1.6H(>i) and by an outer cylinder of radius 
r e = 10 and height — 1.6H(r e ) < z < 1.6H(r e ). These surfaces 
will be indicated as S; and S e respectively. Above the disk the 
control volume is closed at the surface determined by the height 
scale 1.6H(r) between r, = 1 and r e = 10. This surface will be 
indicated as S s . 

The four cases shown in Fig. [8] are characterized by a differ- 
ent outer accretion rate (calculated through S e , i.e. at r e = 10): 
when the poloidal resistivity parameter a m or the anisotropy co- 
efficient Xm are increased the outer accretion rate, indicated as 
M ae , becomes smaller. Moreover it is evident that in the cases 
characterized by a m = 0.1 (dotted and dash-dotted line) the 
accretion rate is at least one order of magnitude higher than 
the initial one (see Table [TJ which was determined to balance 
the advection and the diffusion of poloidal field lines inside the 
disk: it is therefore clear that in these cases the advection of the 
field should dominate thus explaining the peculiar dynamical ef- 
fects of the a m = 0.1 simulations shown in the previous Section. 
Despite of having such a high accretion rate, we recall that the 
low resistivity case at low resolution [dash-dotted line) does not 
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Fig. 8. Radial dependence of the accretion rates at t - 400. 
The plots refer to the cases with the Ohmic heating suppressed 
(/ = 1): (ar m = l,Xm = 3 solid line), (ar m = l,Xm = 1 dotted 
line), (a m = 0.1, Xm = 1 dashed line) and (a m = 0.1, Xm = 1 
, low resolution dot-dashed line line). For the two high diffu- 
sivity cases a power-law approximation (see text) is also plotted 
(long-dashed line). 



show these features, thus supporting the idea that the strong ad- 
vection is balanced by a high numerical diffusion. 

Besides of being characterized by a different M ae , the solu- 
tions show also a different slope of the accretion rate inside the 
"launching region". This slope is clearly linked to the amount of 
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Fig. 9. Temporal evolution of the ejection efficiency. The plots 
refer to the cases (a m — I, Xm — 3 solid line), (a m — I, Xm — 1 
dotted line), (a m — 0.1, Xm — 1 dashed line) and (a m — 0.1, 
Xm — 1 > low resolution dot-dashed line line). 



mass that is extracted from the disk and goes into the outflow: 
defining a simplified power-law radial behavior for the accretion 
rate (Ferreira & Pelletier [T995T l 



M a (r) = M ae - 



we can find, imposing mass conservation inside the accretion 
disk, a simple relation between the ejection efficiency and the 
ejection parameter 



A/ ae 



(31) 



Here the ejection efficiency 2Mj/M ae is defined as the ratio of 
the mass outflow from the disk surface 5 S and the outer accretion 
rate calculated at r e = 10; the factor 2 takes into account the two 
bipolar jets. It is easy to see that the higher the ejection efficiency 
the higher the ejection parameter ^ and therefore the steeper is 
the slope of the accretion rate. 

In Fig.|9]we plot the ejection efficiency 2Mj/M ae as a func- 
tion of time for the four simulations taken into account in this 
Section: it is possible to see that after an initial transient, this 
efficiency reaches an almost constant value, suggesting the at- 
tainment of a stationary final state. On the other hand the solu- 
tions reach different efficiency values, increasing from 20% for 
the (a m - \, x m - 3) case to 55% for the (a m = 0.1, Xm = 1) 
simulation. This result is consistent with the plot shown in Fig.Hl 
where we showed that the forces acting vertically inside the disk 
change sign at lower height scales in less diffusive cases: since 
the disks become denser towards the midplane this determines 
also a higher outflow rate. As usual the low resistivity case at 
low resolution shows an anomalous behavior (dash-dotted line) 
having an ejection efficiency typical of the a m — 1 cases. 

Moreover it is clear the direct correlation between the ejec- 
tion rate and the ratio \B^,\/B p plotted in Fig. [5] This means that 
the main mechanism accelerating the outflow depends strongly 
on the outflow rate: the less mass loaded jets are efficiently cen- 
trifugally accelerated, while in the more loaded jets the higher 
inertia strongly bends the magnetic field in the azimuthal direc- 
tion right above the disk surface; as stated by Eq. ( 130) . the out- 



flow is therefore mainly accelerated by the pressure gradients of 
the toroidal field (see also Anderson et al. 12004) . 

Consistently the solutions characterized by a higher ejection 
efficiency show a steeper slope of the accretion rate as shown 
by Eq. (f3TT >: as an example, we plotted in Fig. [8] the power-laws 
calculated solving Eq. OTb for £ for the cases (a m = l,Xm = 3) 
and (a m = l,Xm = 1) (long-dashed lines). The ejection param- 
eters £ calculated for these two cases are £ = 0.09 and £ = 0.17 
respectively. 

5. Angular momentum transport 

It is commonly accepted that disk-winds are a viable mechanism 
to extract angular momentum from accretion disks without re- 
curring to viscous torques to allow accretion. In this Section we 
show how the bipolar jets can extract angular momentum from 
the underlying disk and through which channel. 

The angular momentum transport is regulated by the angular 
momentum flux through the surface delimiting the control vol- 
ume defined in Section [4] We can therefore define an accretion 



torque «/ acc — >/ a cc,kin 

terms: 



+ A 



defined by the sum of the two 



' acc.kin 



and 



I ripw^H • dS - I , 



r e pu,pu ■ dS 



J at 



f nB^S-dS- \ r e B 



B dS 



(32) 



(33) 



For the accretion torque we use a positive sign if it increases the 
angular momentum contained inside the control volume and a 
negative sign if it extracts it. The term Acc.kin defines the flux 
of angular momentum inside the control volume due to the ac- 
cretion motion, while A C c,mag is the magnetic torque which can 
transport angular momentum along the radial direction inside 
the disk. This magnetic torque represents a small contribution 
to 7 acc : in all of our simulations the ratio Acc,mag/ Action is al- 
ways negligible, around -10%. A small fraction of the disk angu- 
lar momentum is extracted by the magnetic torque also radially, 
thus helping accretion. On the other hand, the main contribution 
to Acc.kin is given by the integral on the outer radius, where both 
the accretion rate and the specific angular momentum of the disk 
are higher. Assuming a Keplerian rotation, a reasonable approx- 
imation for A cc is therefore given by (see also Eq. d28ll): 



Acc ~ M ae yjGMr e 



(34) 



In the same way we can define the torque exerted by the 
outflow on the disk 7j et = /jetjrin + /j e t,mag as the sum of the flux 
of angular momentum 7j e t.kin and of the magnetic torque 7j e t,mag 
at the disk surface: 



jet,kin 



I 



rpu^u ■ dS 



and 



jet, mag 



= f 

ag I 

Js« 



rB.B ■ dS 



(35) 



(36) 



For the jet torque we use a positive sign if extracts angular mo- 
mentum from the disk. 

In a steady situation the accretion and the jet torque should 
be equal: 27j et = 7 acc . Using the approximation made in Eq. ( |34l , 
the relation 



M ae VGMr e ~ 27j 



jet 



(37) 
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Fig. 10. Temporal evolution of the jet torque 2/j et (left panel), of the ratio between the jet and accretion torque 2jj el /J. dcc (central 
panel) and of the ratio between the jet magnetic and mechanical torque /j e t,mag//jet,kin (right panel). The plots refer to the cases 
(a m — l,Xm — 3 solid line), (a m — l,Xm — 1 dotted line), (a m — 0.1, Xm — 1 dashed line) and (a m — 0.1, Xm — 1 , low resolution 
dot-dashed line line). 



shows clearly that the accretion rate at the outer radius of the 
launching region r e is mainly controlled by the torque exerted 
by the outflow on the accretion disk. In the left panel of Fig. [10] 
we plot the temporal evolution of the total torque 2jj et exerted 
by the bipolar jets of the four simulations performed without 
Joule heating: it is clear the correlation between the outer ac- 
cretion rate shown in Fig. [8] and the jet torque plotted here, as 
stated by Eq. d37l ). The high accretion rate of the a m — 0. 1 cases, 
which determines the peculiar behavior of these solutions, is de- 
termined by the high torque exerted by the jet, at least ten times 
higher than the torque needed to maintain the initial rate (see 
TableQ]). 

The curves plotted in the left panel of Fig. [TUl show a slow 
decrease in time, thus suggesting that our solutions did not reach 
a final steady state. On the other hand, the ratio 2/j et / 7 acc , plotted 
in the central panel of Fig.[10j stays approximatively constant in 
time and equal to one, after an initial transient: this shows that 
the accretion rate of the disk reacts quickly to the slowly varying 
jet torque suggesting that our simulations are evolving through a 
series of quasi-stationary states. 

Even if the outflows of the different simulations are exerting 
almost the same torque (left panel, Fig. [Tol l, being a factor two 
smaller just in the case characterized by a high and anisotropic 
magnetic diffusivity (solid line), the ratio between the magnetic 
-/jet.mag and the mechanical torque ~/j e t,kin shows noticeable dif- 
ferences (right panel in Fig. ITOb . The ratio goes from a value 
around unity for the less dissipative case (a m = 0.1, Xm = 1) up 
to a value > 4 for the (a m = l,Xm - 3) simulation. Once again 
we point out that the low resistivity simulation at low resolution 
(dot-dashed line) behaves like a higher resistivity case, show- 
ing /j e t,mag/-/jet,kin ~ 4. This quantity has an easy interpretation: 
the angular momentum extracted from the disk is in fact stored 
initially in the toroidal magnetic field and it is then transferred 
starting from the disk surface to the outflowing plasma which 
is therefore centrifugally accelerated. The ratio /j e t,mag/-/jet,kin is 
therefore a measure of the efficiency of the magneto-centrifugal 
mechanism: the higher is this ratio, the more specific angular 
momentum is available at the disk surface, centrifugally acceler- 
ating the plasma to higher poloidal terminal speeds (see Section 



6. Energy budget of the disk-jet system 



Table 2. Values of the different contributions to the total accre- 
tion power £ acc calculated at t = 400. The simulations are char- 
acterized by their diffusivity parameters a m and^ m ). 



Simulation 


^acc.mec / ^acc 


^act^mag/ ^acc 


^acc^hm / ^acc 


(0.1,1) 


0.97 


0.14 


-0.11 


(0.1,1) low r. 


1.29 


0.09 


-0.38 


(1,1) 


0.97 


0.10 


-0.07 


(1,3) 


0.99 


0.05 


-0.04 



consider the energy fluxes through the surfaces of the control 
volume defined in Section [4] We then define an accretion power 
£acc = £acc,mec +£acc,ma g + £acc,thm given by the sum of the terms: 

£acc,mec = j[ (j + ' dS ~ jT (j + #g)p« ' ^ (38) 

£acc.mag = f E X B ■ dS - f E X B ■ AS (39) 

Js, Js L . 

£acc,thm = f ~^PU AS- { -^~PU ■ dS (40) 

where £acc,mec, £acc,mag and £acc,thm represent the mechanical (ki- 
netic + gravitational), Poynting and enthalpy flux through the 
cylindrical surfaces at the external and internal radius of the 
launching region. For the accretion power we will use a posi- 
tive sign if it increases the energy contained inside the control 
volume. 

In the left panel of Fig.[TT]we plot the temporal evolution of 
the total accretion power while in Table [2] we show the different 
contribution to E^cc at f = 400 for the cases taken into account 
in this Section. We can see that the dominant term is always the 
mechanical power liberated in the accretion: the Poynting flux 
increases slightly the total energy contained in the control vol- 
ume, while the enthalpy flux always acts to decrease the accre- 
tion power, advecting the thermal energy at the inner radius. It 
is important to notice that the simulation with a low magnetic 
diffusivity performed at a lower resolution (second line) shows 
a higher enthalpy flux due to the higher numerical dissipation. 
Neglecting the enthalpy and the magnetic contribution and as- 
suming a Keplerian rotation, the energy liberated by accretion 
can be safely approximated by (see also Eq. (1281): 



In this Section we finally study the energetics of the disk-jet sys- 
tem. In order to study the energy balance of our simulations we 



(41) 
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Fig. 11. Temporal evolution of the accretion power E^cc (left panel); the ratio between the jet and accretion power 2£ , j et /£ acc and 
the ratio between the radiated and accretion power 2£j e t/£acc (central panel); the ratio between the jet magnetic and kinetic power 
£jet,mag/£jet,kin (right panel). The plots refer to the cases (a m — 1, Xm — 3 solid line), (a m — 1, Xm — 1 dotted line), (a m — 0.1, 
Xm — 1 dashed line) and (a m =0.1, Xm — 1 > l° w resolution dot-dashed line line). 

Table 3. Values of the different contributions to the total jet power £j et calculated at t — 400. The simulations are characterized by 
their diffusivity parameters a m and^ m ). 



Simulation 

(amvfm) 


£jet,kin / ^jet 




^jet,mag/ -^jet 


£jet,thm /^jet 


(0.1,1) 


0.42 


-0.73 


1.28 


0.03 


(0.1,1) low r. 


0.13 


-0.24 


1.05 


0.06 


(1,1) 


0.25 


-0.49 


1.22 


0.02 


(1,3) 


0.13 


-0.26 


1.12 


0.01 



)0 



1.00 
S 0.95 
"i 0.90 
'» 0.85 
0.80 
0.14 

„ o.iaf- 

j 0.10 
~i 0.08 
. S 0.06 
" 0.04 
0.02 
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time 
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This approximation shows clearly that the liberated power is 
mainly determined by the inner accretion rate: the values of E- dcc 
of our simulations at t = 400 (left panel of Fig. fTTb are clearly 
correlated with the inner accretion rates visible in Fig. [8] 

To characterize the energy extracted by the outflow we define 
a jet power E jet = E ietMn + £j e t, g rv + E ieUmag + Ej eUhm given by the 
sum of the fluxes of kinetic, gravitational, magnetic and thermal 
energy at the disk surface 5 S : 



^jet.kin — 

^jet,grv — 

F- - 

^jet.mag 

^jet.thm = 



Impute 
j «DgP« dS 

f — 



E x B ■ AS 



Pu ■ AS 



(42) 
(43) 
(44) 
(45) 



For the jet power we will use a positive sign if it extracts energy 
from the control volume. 

Since in this Section we are taking into account cases char- 
acterized by / = 1, we can also define the radiated power E coo \ 
as 



J cool 



f A CO oidV= f f,J 

Jv c Jv c 



JAY 



(46) 



where V c is the control volume. In a steady situation the en- 
ergy fluxes should balance to give 2£j et = £ a cc - ^cooi- In me 
central panel of Fig. [TT] we plot the time evolution of the ratios 
2Ej et /E acc and E coo \/E acc . We can see that most of the accretion 
power is liberated in the jet, increasing from ~ 90% for the more 
dissipative case (a m = 1, Xm = 3) up to ~ 98% for the less 
dissipative one (a m - 0.1, Xm = 1)- Correspondingly the radi- 
ated power decreases from ~ 10% down to ~ 2% respectively. 



Anyway the two contibutions sum up to approximatively one, 
suggesting a quasi-stationary situation. 

As it was noticed for the jet torques, also the jet powers do 
not show a huge difference between the simulations. On the other 
hand the different contributions to the total jet power depend 
more clearly on the simulation parameters. In Table [3] we see 
that the greater contribution to the jet power at the surface of 
the disk comes from the Poynting flux, since most of the en- 
ergy is stored in the toroidal magnetic field. The enthalpy flux 
is always negligible: our jets are cold, undergoing a continuous 
adiabatic expansion starting from the disk surface. The kinetic 
and the gravitational fluxes are strongly correlated with the out- 
flow rates of the different cases: the gravitational contribution is 
negative since at the disk surface the outflow is still inside the 
potential well of the central object. 

In the right panel of Fig. [TT] we plot the temporal evo- 
lution of the ratio between the jet Poynting and kinetic flux 
£jet,mag/£jet,kin- A higher £j et ,mag /£jet,kin ratio indicates that more 
magnetic energy per particle is available at the disk surface to 
accelerate the outflow. The Poynting flux available at the disk 
surface is then converted into kinetic energy along the outflow 
and at the upper end of the computational box the kinetic flux 
is dominant, as confirmed by the fact that the flow is super-fast- 
magnetosonic. The asymptotic speed of the outflow is therefore 
higher for cases characterized by a higher £j e t,mag/£jet,km, as con- 
firmed by Fig. [12] where we plot an average value of the poloidal 
speed normalized to the toroidal (Keplerian) speed at the base 
of each streamline. The poloidal velocity measured at the upper 
end of the computational domain is a few times the Keplerian 
speed at the base of the streamlines, consistently with the fact 
that different classes of astrophysical jets show velocities of the 
same order of the escape velocity from the central object (Livio 
fT999l . 
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Fig. 12. Average poloidal speed u po i along the jets calculated at 
t = 400 normalized to the value of the toroidal (Keplerian) speed 
u^o at the base of each streamline. These curves have been ob- 
tained averaging the ratio u po i/u^o on all the streamlines outflow- 
ing from the upper boundary of the computational box. The line 
style used refers to the same cases as in Fig.fTTI 
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Fig. 14. Rotation rate of the magnetic field calculated along the 
fieldlines whose footpoints are located at r = 2 {solid line), 4 
{dashed line) and 8 {dot-dashed line. 



will analyze the simulation characterized by a m = \,Xm — 3 and 
/ = 1 : besides of showing some features of a stationary solution, 
its parameters a m ,Xm> I 1 an d e are typical of the cold self-similar 
steady solutions found by Casse & Ferreira (2000a), allowing a 
direct comparison between steady and time-dependent solutions. 

We recall a few relations which are valid for a steady ax- 
isymmetric solution of the ideal MHD equations. The poloidal 
speed Hp is related to the poloidal magnetic field B p 



K{W) 
Up = B p 

P 



(47) 



where K ( V F) is a constant along every field line marked by the 
flux function *F and it represents the ratio between the mass flux 
and the magnetic flux. This relation clearly states that in a steady 
situation the poloidal speed and magnetic field are parallel. The 
toroidal speed 11$ is related to the rotation rate of the magnetic 
surface £2(*P) by 



iu=OC¥)r+ ^-B, 



(48) 



Fig. 13. Angle between poloidal speed and poloidal magnetic 
field for the case characterized by a m = l,Xm = 3 and with the 
Ohmic heating suppressed. This quantity is calculated at t — 400 
along three sample field lines whose footpoints are located at ro 
= 2 {solid line), 4 {dashed line) and 8 {dot-dashed line). 



7. A quasi-steady solution 

In this Section we will try to characterize one of the simula- 
tions performed as a steady solution. We already pointed out 
that none of our solutions are perfectly steady: the low diffu- 
sivity case shows a continuously evolving magnetic structure, 
while all the simulations show a slowly evolving jet torque (left 
panel in Fig. |T0T > and accretion power (left panel in Fig. ITTb . On 
the other hand the constant values assumed by the accretion ef- 
ficiency (Fig. |9), by the ratio between the jet and the accretion 
torque (central panel in Fig. |T0T > and by the ratio between the 
jet and the accretion power (central panel in Fig. [TTT i. indicate 
that our solutions, at least those with or m = 1, are slowly evolv- 
ing through a series of quasi-stationary states. In this Section we 



The total specific angular momentum IQ¥), which is constant 
along each field line, is given by 



Iffl = nu - 



rB 4 
KQV) 



(49) 



where r& is the Alfven radius. Each field line is therefore char- 
acterized by the non-dimensional constants k, or mass loading 
parameter, and A, which gives a measure of the magnetic lever 
arm r A /r : 



00 



A 



Qr 



(50) 
(51) 



where Tq is the cylindrical radius of the footpoint of the field line 
and Bq is the value of the magnetic field at rn. 

We therefore tried to characterize our simulation computing 
these steady invariants along three sample field lines whose foot- 
points are located at ro = 2 {solid line in the next Figures), 4 
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Fig. 15. Mass loading parameter calculated along the fieldlines Fig. 16. Magnetic lever arm parameter calculated along the field- 
whose footpoints are located at r = 2 (solid line), 4 (dashed lines whose footpoints are located at r$ = 2 (solid line), 4 (dashed 
line) and 8 (dot-dashed line). line) and 8 (dot-dashed line). 



(dashed line) and 8 (dot-dashed line). In Fig. [13] is plotted the 
angle between the poloidal field and the poloidal speed: while 
inside the disk, where the magnetic resistivity is effective, the 
accreting plasma flows perpendicularly to the field lines, in the 
outflow, where the ideal MHD holds, the matter flows along the 
lines. In Fig. [14] is plotted the Q. invariant: the rotation rate of 
the field lines, defined by a reference frame in which the electric 
field is zero, is close to the Keplerian rate, since the field lines 
are anchored in the accretion disk and corotate with it. Moreover 
it is possible to see that Q is slightly higher at the base of the out- 
flow: as in the low resistivity case but without noticeable conse- 
quences, the footpoint of the field line is advected towards the 
central part of the disk, rotating a little bit faster than its outer 
end. 

The mass loading parameter k is plotted in Fig[l5] while the 
A invariant is shown in Fig. [16] The two quantities show a much 
more constant behavior along the inner field line (solid line), 
indicating that the central part of the outflow has reached a more 
stationary state, while the outer part of it is still slowly evolving. 
Anyway a clear trend emerges in these plots: the mass loading 
decreases going from the inner to the outer part of the outflow 
while the magnetic lever arm increases. The relation between 
the k and A parameters is shown in Fig. [T7] where we plotted 
the values of these quantities calculated at the Alfven point on 
different fieldlines anchored in the disk. The higher values of k 
are found on the inner (r ~ 1) fieldlines while the lower one in 
the outer part (r ~ 10) of the launching region. For comparison 
it is also plotted the relation: 



10F 



(52) 



found by Weber & Davis ( 1967) for a radial wind geometry. The 
rough behavior A oc k~ 2 ^ is found also by Ouyed & Pudritz 
( 1997) and Anderson et al. d2004l >. 



We point out that in cold (adiabatic) steady solutions the fol- 
lowing relation is valid (see e.g. Blandford & Payne 1982) along 
each field line: 



^jet.mag 



^jet,kin 



2^=2(M) 



(53) 
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Fig. 17. k-A relation (squares). The values are calculated at the 
Alfven point along different fieldlines anchored in the launching 
region of the disk (1 < r < 10). The value of k increases going 
from the outer to the inner fieldlines. For comparison it is also 
plotted the relation Eq. d52l found by Weber & Davis ( 1967} for 
a radial wind geometry (dot-dashed line). 

We plotted an averaged value over the disk surface of these two 
ratios in Fig. [TO] and QT] The magnetic lever arm is therefore a 
measure of the efficiency of the magneto-centrifugal mechanism: 
assuming that all the Poynting flux available at the disk surface is 
completely converted into poloidal kinetic energy, the following 
relation for the poloidal asymptotic speed of the outflow u Pi00 
holds: 



Qr V21 - 3 



(54) 



Therefore the poloidal speed, expressed in units of the rotation 
speed at the footpoint, asimptotically assumes higher values on 
the outer field lines. 

Finally, it is interesting to compare our solution with the cold 
self-similar steady solutions found by Casse & Ferreira (2000a) 
characterized by the same parameters a m , Xm, £ and fi. The simu- 
lation presented in this section shows greater values of the mass 
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loading k and of the ejection parameter £ ~ 0.1 (calculated in 
Section |H and lower magnetic lever arms. The steady solution, 
characterized approximatively by k ~ 2 x 10~ 2 , £ ~ 10~ 2 and 
A ~ 35, has therefore a much lower outflow rate and higher ter- 
minal speed. 

A possible physical explanation is given by the shape of 
the current circuits shown in Section [3~2l Ferreira (1997J has 
shown that high ejection efficiencies > 0.5) are obtained if the 
poloidal current enters the disk at its surface and that no steady 
solution crossing all the critical points is possible. Its solutions 
are always characterized by a lower ejection parameter ^ < 0.5 
and by a current outflowing from the disk surface. In our so- 
lutions the current returns and enters the disk in small region 
1 < r < 2 near the central sink region: due to the internal bound- 
ary itself the rotation of the outflow goes suddenly to zero at the 
inner edge of the jet, creating a current sheet which is forced to 
flow back inside the disk. This is therefore a numerical effect and 
a proper treatment of the region of interaction between the disk 
and the central object is required to understand how the currents 
close in the inner region of the system. Anyway it is possible 
to speculate that the returning current increases the ejection ef- 
ficiency as proposed by Ferreira ( 1997 ): this could be also con- 
firmed by the fact that in our solution the mass loading parameter 
k increases towards the center, where the current flows back into 
the disk. 

Another possible explanation is determined by the numeri- 
cal dissipation which can increase the thermal energy and the 
pressure in the launching region, increasing the outflow rate as 
shown in Section [3~T1 (see also Casse & Ferreira 2000b). If on 
one hand we showed that a low resolution yields the effects of a 
high disk diffusivity, on the other hand we cannot be confident 
that the standard resolution used in the simulations is sufficient 
to accurately resolve the vertical structure of the accretion disk: 
it is well known that it is difficult for the type of algorithms as the 
one used in this paper to handle low densities, as in the launch- 
ing region, where the disk expands a lot. This is also suggested 
by the behavior of the invariants k and A, e.g. on the inner field 
line considered: they show a slow transition to their more or less 
constant values, indicating that the transition between the resis- 
tive and the ideal MHD regimes happens on a scale larger than 
the one determined by the vertical profile of the diffusivity (Eq. 

8. Effects of Ohmic heating 

We finally make a few considerations about the two simula- 
tions performed in which all the dissipated magnetic energy is 
released locally inside the disk (f = 0) and we will make a 
comparison with the corresponding cases in which the dissipated 
energy is radiated away (/ = 1). We recall that these two sim- 
ulations are characterized by (a m = 0.1, Xm = 1) and (a m = 1, 

Xm = 3). 

The dissipative Ohmic term A = i]J ■ J acts now to increase 
the thermal energy and the entropy of the plasma inside the disk. 
This has two important consequences: the energy dissipation at 
the midplane changes the disk thickness and modifies its struc- 
ture while an increase of the pressure gradient at the disk surface 
allows more matter to be loaded in the outflow. 

In Fig.[T8]we plot the time evolution of the ejection efficiency 
of the two "hot" solutions and of the two corresponding "cold" 
simulations: both the "hot" simulations are characterized by a 
higher efficiency; moreover, the case with a m — 1 (solid line) 
shows at the end of the simulation a steep increase. This is due 
both to a decrease of the outer accretion rate and to an increase 
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Fig. 18. Temporal evolution of the ejection efficiency for the sim- 
ulations in which the magnetic energy is dissipated locally as 
Joule heating: (a m = 1, Xm = 3, solid thick line) and a m = 0.1, 
Xm = 1, dashed thick line). With a thinner line we also plot the 
curves which refer to the corresponding cases in which the dis- 
sipated energy is radiated. 

of the outflow rate: the energy dissipation steepens the radial 
gradient of the thermal pressure inside the disk, thus slowing 
down the outer accretion; at the inner radius the gradient is so 
steep that the disk does not accrete anymore and all the matter is 
pushed by the pressure gradient in a highly unsteady outflow. 

The disk thermal energy of the "cold" cases is almost con- 
stant in time. The "hot" case with a m = 0.1 shows a similar 
behavior, with a disk thermal energy just slightly higher than 
the corresponding "cold" simulation. It is likely that the higher 
outflow rate can balance the small increase of the disk thermal 
energy which is equivalent to the 2% of the accretion power 
that was "radiated" in the corresponding "cold" simulation (see 
Section |6j. On the other hand the disk thermal energy of the 
"hot" counterpart of the a m = 1 simulation continuously in- 
creases, not balanced by the energy extracted by the outflow, 
leading to the unsteady behavior shown in Fig. [18] It is impor- 
tant to notice that the behavior of this simulation depends a lot 
on the expression that we used for the magnetic diffusivity (Eq. 
( fT3l ) which is proportional to the disk thermal height scale: the 
Ohmic heating determines in fact an increase of the disk thick- 
ness that in turn implies a higher magnetic diffusivity. 

9. Summary and conclusions 

Our calculations have followed the long-term evolution of an 
axisymmetric quasi-Keplerian magnetized disk up to the estab- 
lishment of an inflow/outflow configuration. The accretion flow 
is driven by extraction of angular momentum of the disk by the 
jet; this is shown by reaching balance between the accretion and 
jet torques. The magnetic torque of the jet is most efficient close 
to the surface of the disk extracting angular momentum from the 
accretion flow; it stores angular momentum in the toroidal mag- 
netic field that then accelerates the outflowing plasma. Therefore 
the simulations have succeeded in demonstrating that the mag- 
netocentrifugal mechanism originally proposed by Blandford & 
Payne can launch jets, provided certain physical conditions on 
the magnetic resistivity and initial field configuration are satis- 
fied. 
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In particular we have shown that an isotropic (x m = 1) resis- 
tive configuration with a m = 0.1, due to the stronger advection 
of the field compared to its diffusion, produces highly unsteady 
magnetic structures, like it was displayed in Fig.|2]and in the left 
panel of Fig. [7] This is in agreement with the stationary models 
of Casse & Ferreira (2000a ) according to which it is not possible 
to obtain a steady outflow with such a low a m parameter except 
for highly anisotropic configurations (x m > 10 2 ). This would 
be in agreement with the results shown in Section [5] a much 
stronger toroidal resistivity should reduce the torque exerted by 
the jet on the disk and thereby its accretion rate (Eq. (l37l>). On 
the other hand we have also shown how this unsteady behav- 
ior depends strongly on the resolution assumed and therefore on 
the numerical dissipation: the same case a m = 0.1 performed 
with a resolution four times lower presents many of the char- 
acteristics of a solution with a m = 1. It can be presumed that 
the quasi-stationary behavior found by Casse & Keppens (2002 
2004) with an isotropic a m = 0.1 parameter can be affected by 
the resolution used. 

On the other hand, also the cases characterized by a m — 1 are 
not perfectly stationary, despite of showing an ordered magnetic 
configuration favorable for a steady launching: in all the simu- 
lations performed some integrated quantities, like the jet torque 
(Fig. [TO]) or the accretion energy (Fig.fTTI). are still slowly evolv- 
ing. Even the case characterized by an anisotropic (x m = 3) dif- 
fusivity, needed by the stationary models, does not show a per- 
fectly stationary behavior; moreover it can be characterized by 
adimensional quantities like k or A which differ a lot from 
the analogous found in the cold solutions of Casse & Ferreira 
(2000a). In Section[7jwe proposed some physical and numerical 
reasons to explain this behavior. 

Nevertheless the constant value assumed in time by some 
quantities, like the ejection efficiency (see Fig. [9]) or the ra- 
tio between the accretion power and the jet energy flux (cen- 
tral panel of Fig. ITTb suggests that our solutions are slowly 
evolving through a series of quasi-stationary states. Some clear 
trend emerges: increasing the poloidal and/or the toroidal resis- 
tivity the ejection efficiency decreases, going from 55% in the 
(a m = 0.1, Xm = 1) case t° 20% for the (a m = l,Xm = 3) sim- 
ulation. From the energetic point of view, our simulations show 
that more than 90% of the energy liberated in the accretion in- 
flow is released in the jet and that less dissipative cases produce 
slightly more powerful jets. Correspondingly, the radiated power 
that must take care of the Joule heating dissipation is less than 
10%. In Section|8]we have also shown that in the more dissipa- 
tive case with a m = 1, if this energy is released inside the disk 
instead of being radiated, the accretion flow is disrupted and a 
highly unsteady outflow is formed. 

We have also demonstrated how the efficiency of the 
magneto-centrifugal mechanism, as measured by the ratio of the 
Poynting flux to the kinetic flux at the disk surface, is affected 
by the resistive configuration. A higher value of the poloidal 
and/or toroidal resistivity determines a greater magnetic lever 
arm, linked to the energy flux ratio as stated by Eq. ((53). Our 
simulations show also that the A parameter, as measured by Eq. 
( f53l ). is related with a good approximation to the ejection param- 
eter £ (in the cases where this parameter has been measured): 

A ~ 1 + ^ (55) 

as stated by the self-similar models by Ferreira (1997). 

We finally try to test if observations of outflows from classi- 
cal T-Tauri stars can give some constraint on the parameters of 
our simulations. The size of the launching region measured in 



our simulations (from 0.1 to 1 AU) is consistent with the es- 
timates derived from observations (Bacciotti et al. 2000). On 
the other hand, the one-sided ejection efficiency inferred from 
the observations gives values, even with high uncertainty, in the 
range 0.01-1 (Cabrit 2002): this result seems to favor our solu- 
tion characterized by (a m = 1, Xm = 3). Even if, in agreement 
with observations, all our solutions are characterized by a speed 
at the upper boundary of the computational domain which is a 
few times the escape speed from the potential well of the cen- 
tral object (Fig.fTZT). a quantitative comparison with the poloidal 
speed (Bacciotti et al. 2000) and rotation signatures (Bacciotti et 
al. l2052l Coffey et al. 120041 Woitas et al. EUD31 > of the observed 
outflows is much more difficult: the size of our computational 
domain, which reaches a physical scale along z equal to 12 AU, 
is around three times smaller than the distance from the source 
investigated by current HST observations (> 30 AU). Anyway, 
as it has been pointed out by Ferreira et al. (120061 ). a solution 
characterized by a magnetic lever arm A ~ 10 successfully re- 
produces the values of both the poloidal and toroidal speeds at 
the currently observed spatial scale: the highest lever arm value 
found in our simulations (A ~ 9) is observed in the outer part of 
the launching region of the case (a m = l,Xm = 3) (Fig.fTTI). 

Two specific conclusions and difficulties of our calculations 
must be mentioned. First, the magnetocentrifugal process ap- 
pears to operate only with relatively strong magnetic fields in 
the disk; a strong field tends to inhibit the development of turbu- 
lence that could be the origin of magnetic resistivity. Instabilities 
(shear, Kelvin-Helmholtz) other than magneto-rotational might 
produce the required turbulence. On the other hand, the resistiv- 
ity inside a protostellar disk can be extremely high between 0. 1 
and 10 AU, where the grains become the dominant charge car- 
riers (Wardle 1997). The coupling between the accreting plasma 
and the magnetic field is therefore strongly, perhaps too much, 
reduced. The second point is that outflow rates measured in our 
simulations turn out to be rather high, which applies well to jets 
in star formation regions, where extended bipolar structures are 
observed; for AGN nuclei one must go to the relativistic limit 
(which we have not done here), and also discuss how the (pos- 
sibly subsonic) wind from the disk outer part interacts with the 
surrounding galaxy. 
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